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CN ■ ABSTRACT 

> ; 

C2 Context. Realistic atmospheric models that link the properties and the physical conditions of supernova ejecta to observable spectra are required 
for the quantitative interpretation of observational data of type la supernovae (SN la) and the assessment of the physical merits of theoretical 
supernova explosion models. The numerical treatment of the radiation transport - yielding the synthetic spectra - in models of SN la ejecta 
in early phases is usually carried out in analogy to atmospheric models of 'normal' hot stars. Applying this analogy indiscriminately leads to 
inconsistencies in SN la models because a diffusive lower boundary, while justified for hot stars, is invalid for hydrogen and helium-deficient 
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f ^ i supernova ejecta. In type la supernovae the radiation field does not thermalize even at large depths, and large optical depths are not reached at 
*»«^. , all wavelengths. 

' Aims. We aim to derive an improved description of the lower boundary that allows a more consistent solution of the radiation transfer in SN la 
and therefore yields more realistic synthetic spectra. 
Q ■ Methods. We analyze the conditions that lead to a breakdown of the conventional diffusion approximation as the lower boundary in SN la. 
!— I For the radiative transfer, we use a full non-LTE code originally developed for radiatively driven winds of hot stars, with adaptations for 
CX) , the physical conditions in SN la. In addition to a well-tested treatment of the underlying microphysical processes, this code allows a direct 

comparison of the results for SN la and hot stars. 
^ Results. We develop a semi-analytical description that allows us to overcome some of the limiting assumptions in the conventional treatment 
i— I . of the lower boundary in SN la radiative transfer models. We achieve good agreement in a comparison between the synthetic spectrum of our 
^ ' test model and an observed spectrum. 
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1. Introduction 

Type la supernovae (SN la) have become an invalua ble tool for 
the determination of the cosmological paramete rs (Riess et al. 
1998[ IPerlmutter etaD Il999t iRiess et all Eooil iTonrv et alJ 
2003) as their exceptional brightness makes them observable 
even at large cosmological distances. Using SN la for dis- 
tance determination requires knowledge of their absolute lu- 
minosities; however, SN la are not perfect "standard can- 
dles" in this respect because they show an intrinsic scatter in 
their properties, in particular in the peak brightness. The ap- 
plication of SN la for cosmology therefore relies on empir- 
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ical relationships between the pe ak bright ness and other ob- 
served characteristics (e.g., Phi mpslll993t lHamuv et alj fl 996 ; 
IRiess et alJ Il996t IPerlmutter etalJll997l) . To first order the 
shape of the light curves and certain spectral properties is 
determined by the m ass of s ynthesiz ed 56 Ni and its distribu- 
tion within the ejecta JNueent et alJl995btfHoflich et alJl995t 
IPinto & EastmanlbOOOh . However, the details of the physical 

processes that cause the observed variation are still unclear 

■ II i 

(see Hillebrandt & Niemever 2000 for a review). This uncer- 
tainty is an essential problem because the different calibration 
methods used to derive the peak brightness yield partly dif- 
ferent re^ultSjJmnlyjrig an unaccounted source of systematic 
error (Leibundaut 2004). In addition, the application of SN la 
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for cosmological distance measurement relies on the crucial as- 
sumption that objects observed at high redshifts have the same 
properties as objects in the nearby universe, so that the same 
calibration method for the luminosity differences can be ap- 
plied to all objects. It is, therefore, of fundamental interest to 
derive a physical model that can explain the explosion mecha- 
nism of SN la in detail, including the observed intrinsic vari- 
ability. This will allow a more reliable estimate of the system- 
atic errors of the distance measurement. 

The currently favored models for the explosion mechanism 
of SN la involve the thermonuclear explosion of a carbon- 
oxygen white dwarf (WD) in a binary system. Two progeni- 
tor scenarios are generally considered: the "single degenerate" 
scenario and the "double degenerate" scenario. In the single 
degenerate scenario the WD accretes mass from a red-giant 
companion star. When the WD reaches a mass close to the 
Chandrasekhar mass (Men ~ 1 AM©) the compressional heat- 
ing at the center of the star triggers carbon burning. After a pe- 
riod of a few thousand years of quiet burning a thermonuclear 
runaway disrupts the star Clben & Tutukov 1984; Webbink 
ll984llWooslev et alJ20Q4tlHan & PodsiadlowskJ20Q4h . In the 
"double degenerate" scenario the companion star is also a WD 
with a mass such that the total mass of the system exceeds Mch- 
Due to energy loss by gravitational waves the orbital separation 
of the binary system gradually decreases, leading to a merger, 
whic h triggers the thermonuclear runaway that explodes the 
star dWhelan & Ibenlll973l lNomotol ll982). Potential progen- 
itor systems have been found in recent years; their numbers, 
however, are too low to explain the obse rved rates of SN la 
JCappellaro et alJl999h.fPauldrachl|2005h suggests on basis of 
the results of Paul drach et al.l Il2004h a potential connection of 
SN la progenitors to a subgroup of central stars of planetary 
nebulae (CSPN). 

The details of the explosion process itself are also still 
subject to a lively debate. The general picture is that a sub- 
sonic deflagration wave ("flame") is ignited near the center of 
the star. The flame travels outward, burning part of the star 
to nuclear statistical equilibrium (NSE). Because the flame 
propagates subsonically, the star can expand while undergo- 
ing burning. This allows partial burning of C and O to in- 
termediate mass elements (Si, S, Mg, Ca) which dominate 
the composition. In contrast to the deflagration, a prompt, 
supersonic detonation of the star would generate predomi- 
nantly iron group elements. This outcome contradicts the ob- 
served composition. No agreement has been reached in the 



the end of the explosion ( Nomoto et al 


1984; 


Wooslev et al. 


1984 iNiemever & Hillebrandt 1995: Reinecke et alJ 12002: 


Ropke & Hillebrandt 2005; Ropke 2005 


or if a (yet unknown) 



mechanism triggers the deflagration to turn into a supersonic 
detonation toward t he end of t he explosi on (delayed deton ation 
transi tion, DDT) faoflich & Khokhlovl fl996t llwamoto et all 
119991) . At present, the resulting composition of DDT models 
seem to favor the latter scenario; however, those models are 
not self-consistent and depend on ad hoc assumptions of the 
occurrence of the DDT. 

Judging the validity of numerical explosion models re- 
quires a comparison of the observational consequences pre- 



dicted by the explosion models to the de facto observations of 
SN la. Realistic radiative transfer models provide the crucial 
link between explosion models and observations. It is possible 
to predict the observational implications of the hydrodynami- 
cal models only if the radiative transfer models are sufficiently 
realistic and detailed. Conversely, such radiative transfer mod- 
els make it possible to establish constraints on the composition 
and structure of SN la from the spectroscopic interpretation of 
observed spectra. 

SN la in phases before and shortly after maximum exhibit a 
spectrum dominated by a few very broad absorption features in 
a non-thermal continuum. These absorption features are mostly 
due to blends of several lines, while the "continuum" itself is 
formed by the overlap of a large number of Doppler-broadened 
metal lines. The true continuum opacities and emissivities that 
determine the overall shape of the spectrum in stars are of mi- 
nor significance in supernovae. At later epochs, the "pseudo- 
continuum" photosphere recedes deeper and deeper into the 
ejecta and eventually disappears as the ejecta become trans- 
parent. Unlike stars, the radiation in SN la is generated within 
the expanding medium itself by the deposition of the energy of 
r-photons that resu lt from the decay chain 56 Ni— » 56 Co— » 56 Fe 
Colgate & McKeell969l) . 

The early epochs where a photosphere is still present are re- 
ferred to as the photospheric phase. In the photospheric phase 
the ejecta can be treated in analogy to hot stars as expanding, 
extended atmospheres. Only the radiative transfer in the photo- 
sphere and in the outer envelope above the photosphere needs 
to be considered. In this setup a steady state is assumed because 
the photon escape timescale in the thin medium generally will 
be much shorter than the expansion timescale. 

In principle, a complete radiative transfer model for SN la 
would require consistent, time-dependent solutions of the pop- 
ulations of all atomic levels and the continuum and line trans- 
fer, including the treatment of energy deposition by the decay 
products of 56 Ni and 56 Co. 

Current models implement various simplifications accord- 
ing to the specific model's purpose. These simplifications 
are necessary because the solution of time-dependent radia- 
tive transfer in three dimensions, including the full coupling 
of radiation and matter, is not yet feasible and some of the 
terms involved in such a consistent solution are shown to be 
or regarded to be of second order. In recent decades syn- 
thetic spectra of SN la have been modeled by several groups 
with a variety of approaches involving different leve ls of 
complexity depending on the application dBranch etal.lf l985: 
i Mazzaii et alJll993| llVIazzali & Lucvlll993t lEastman & Pinto l 
|1 993 f iHoflich et alJl 995UNugent et al.ll 995a|IPaiildrach et al J 
1996; Nugen t et alJll997t Lentz et alJ 1200 ll iHoflichl 120051 
Istehle et all2005UBaron et alJ2006HKasen et all2006l) . 

Highly parametrized models, which implement a simplified 
treatment of physical processes to achieve short run-times seem 
to be suitable for the comparative analysis of a large number 
of observed spectra, while more realistic models are required 
for a deeper understanding of the physical effects leading to 
specific observed properties. In particular, judging the validity 
of hydrodynamic explosion models, as mentioned above, can 
only be performed using radiative transfer models that include 
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a very detailed treatment of relevant physical processes. Such 
detailed models may also be used to validate or invalidate spe- 
cific simplifying assumptions used in less elaborate models. 

In spite of the analogy mentioned before, SN la ejecta differ 
from ordinary stellar atmospheres in several important aspects. 
Techniques and methods that are adopted for stellar atmosphere 
modeling must be carefully checked to verify whether the ap- 
plied approximations are still justified for SN la. In this work 
we present an improved description for the lower boundary of 
the radiative transfer calculation. This is required because the 
assumption of thermalization at large depths, justified for stel- 
lar atmospheres, breaks down for SN la. 

In the work we present here, we us e the computer program 
WMbasic JPauldrach et aljfl994l 1200 lh to obtain a consistent 
solution of the full non-LTE rate equations and a detailed ob- 
server's frame solution of the radiative transfer. This code was 
originally designed for the analysis of the spectra of hot stars 
with radiativel y driven winds, but an earlier version has already 
been used bv lPauldrach et aL dl996h (Paper I) to quantitatively 
investigate the effects of line blocking in SN la. While they 
used a consistent treatment of line blocking, the back-reaction 
of the line opacities on the temperature structure (line blanket- 
ing) was not taken into account. In our present work the current 
version of the code has been further adapted to treat the radia- 
tive transfer in supernovae in a more sophisticated way. 

Section [2] describes the setup of the non-LTE model and 
introduces the numerical scheme used to solve the radiative 
transfer. In Section [3] the physical conditions in the "pseudo- 
photosphere" of SN la are discussed with respect to the solu- 
tion of the radiative transfer and are compared to the situation 
in normal stars (i.e., stars that have a well defined photosphere). 
Section|4]describes the derivation of an improved treatment of 
the inner boundary for the numerical solution of the radiative 
transfer. The results are discussed in Section[5] and a compar- 
ison of a model spectrum with an observed SN la spectrum is 
shown. The conclusions are provided in Section|6] 

2. The radiative transfer model 

The code WMbasic has been successfully used to model ex- 
tended, radiation-driven stellar atmospheres, assuming a ho- 
mogeneous, stationary, spherically symmetric outflow. Here we 
provide a brief outline of the main concepts; details relevant to 
the derivation in Section 0] will be discussed as well. A more 
compreh ensive description of t he numerical m ethods can be 
found in Pauld rach et ail feOOlb and lPauldrachl J2003I) and ref- 
erences therein. 



2.1. General setup 

To derive synthetic spectra of supernovae in early phases the 
analogy to hot stars with extended atmospheres can be used to 
apply similar concepts for the solution. The radiative transfer 
model for supernovae requires the following input: 

- Hydrodynamic structure — This defines the relationship 
between radius and density and velocity at a given epoch 
(i.e., time after explosion). Because the expansion is ballis- 



tic, a homologous structure of the ejecta is generally a safe 
assumption (to describe properties of the ejecta, the veloc- 
ity is usually used as a radial coordinate because it is inde- 
pendent of the epoch). Thus, radius r, epoch t, and velocity 
v are related via r = vt, and a density po(r, to) given at an 
epoch fo scales with t as p(r, t) = po((to/f)r,to)(to/f) 3 , The 
relationship between density and radius also sets the op- 
tical depth scale, which defines the "photospheric radius." 
The hydrodynamic quantities r, v, and p remain fixed dur- 
ing the calculation of the non-LTE model. 

- Luminosity — In SN la the luminosity results primarily 
from the energy that is deposited by y-photons originating 
from the decay of radionuclides synthesized in the explo- 
sion. Therefore, part of the luminosity will depend on the 
radial distribution of the respective elements in the ejecta. 
At the inner boundary of the computational volume an in- 
coming luminosity must be specified to account for the 
radiative energy deposited below that boundary. This lu- 
minosity represents both the instantaneously released pho- 
tons as well as the photons that have been generated at 
earlie r epochs and trapped by the large opacities. (See, 
e.gJArnetllll982tlHoflich & poldilovll99alNugent et alJ 
1997; Pinto & Eastman 2000.) At the current stage of our 
project, the explicit energy deposition within the ejecta has 
not been taken into account; instead, the total luminosity is 
assumed to originate from below the boundary of the com- 
putational grid. Also, the luminosity is not derived from the 
mass of 56 Ni predicted by a specific explosion model but is 
considered as an independent input parameter. 

- Epoch — For homologously expanding ejecta, the epoch 
determines the absolute radius and density scale. Because 
the exact time of explosion of an observed supernova is 
generally not known, the epoch represents an important fit 
parameter. 

- Composition — The elemental abundances can be either 
taken from the nucleosynthesis calculations (to study the 
properties of a specific explosion model) or adjusted inde- 
pendently to fit an observed spectrum in order to determine 
the abundances of the object. In general, the composition in 
supernova ejecta is a function of radius. 

The computation of a self-consistent non-LTE-model re- 
quires the simultaneous solution of a number of physical equa- 
tions, in particular the radiative transfer, the rate equations, and 
the energy equation. Fig. Ogives an overview of the physical 
equations and the physical quantities by which these equations 
are interconnected. 

The general concept behind our code, WMbasic, is to first 
obtain a rough solution with a fast, approximate method, and 
then, based on this solution, obtain a completely consistent so- 
lution with an exact, detailed method. The approximate method 
should provide a solution that is sufficiently close to the final 
solution so that only a few iterations with the much more time- 
consuming detailed method are necessary. 

Our fast, approximate method is based on a Doppler- 
broadened sampling technique for line opacities and emissivi- 
ties. The idea behind this method is to solve the radiative trans- 
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Hydrodynamics 

(from explosion model) 
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Fig. 1. Overview of the physical equations that must be consistently solved for the non-LTE model. The input is fixed by an 
external explosion model. In the middle-left box the rate equations which determine the occupation numbers «, of the atomic 
levels are given. The middle-right box shows the radiative transfer equation for the radiation field. The lower-right box gives 
the energy equation that fixes the temperature structure within the atmosphere. The energy deposition by y photons, shown in 
the lower-left box, is currently taken into account only approximately. This means that the models rely on the input of the total 
luminosity L at the lower boundary of the computational grid instead of the 56 Ni-mass. 



fer for a representative sample of frequency points in the rele- 
vant spectral range. 

In the final iterations, a detailed radiative transfer method 
that does not suffer from the approximations of the first itera- 
tion cycle is used. This method uses an exact observer's frame 
solution, equivalent to a comoving frame solution, which cor- 
rectly treats the angular variation of line opacities and emissiv- 
ities. The line profiles are spatially resolved. Multi-line effects 
and back-reactions of the line opacities on the model structures 
are treated correctly. 

The temperature structure is in practice obtained from bal- 
ancing energy gains and losses to the electron gas (heating 
and cooling rates). This description is equivalent to the con- 
dition of radiative equilibrium (indicated in the energy equa- 
tion in Fig. 0, but is numerically advantageous for physical 
conditions where the opacity is dominated by scattering events 
that do not couple th e radiation field to the thermal pool (see 



mat do not couple tne 
|PauldrachetalJ200lh . 



2.2. Solution of the radiative transfer 

To introduce the nomenclature and equations used later, the so- 
lution applied to solve the radiative transfer in the observer's 
frame employing the fast opacity sampling method is recapitu- 
lated in this section. 

To determine the radiation field that enters into the 
Thomson emissivity, an iteration alternating between the ray- 
by-ray solution and the angular-integrated moments equation 
is performed. Both systems are solved with a Feautrier-type 



scheme (Feautrier 1964) as discussed, e.g., injjviihalas ill 9781) . 
For each frequency point, the iteration is performed twice: first 
for a pure continuum model and afterwards for the full prob- 
lem with continuum and lines. The solution is carried out in 
the usual Cartesian-like p-z-coordinate system where each p- 
ray at a given radius shell corresponds to a p — cos ^-direction 
in spherical coordinates (see Fig.|2j. The transfer equation for 




Fig. 2. p-z coordinate system used to solve the transfer equation 
in spherical symmetry. R denotes the radius of the inner core. 



each p-ray is rewritten for the intensities in positive and nega- 
tive z-direction 



^ = ±(S V -/ V ± ) with dT v = - Xv dz. 
dr v 



(1) 
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These equations have to be evaluated at each frequency point 
of the computational grid. In the remainder of this section the 
index v is omitted for brevity. By introducing the new vari- 
ables u — i (7 + + 7~) and v = 4 (7 + - 7~) these two first order 
differential equations with a single boundary condition (initial 
value problems) can be converted into a second order differen- 
tial equation with a boundary condition for each side: 



Au 
d7 



dv 
d7 



A 2 u 
At 2 



(2) 



where S is the source function. 

To close the system, suitable boundary conditions must be 
specified. At the outer boundary the condition 7 =0 (no radi- 
ation incident from outside) leads to 

Au 



u — v 



At 



— u. 



(3) 



At the inner boundary one has to distinguish between those p- 
rays that intersect the photospheric core (p < R) and those that 
do not (p > R). (See Fig.|2]) For core rays the incident intensity 
has to be explicitly specified, 7 + = 7 C0le , while for non-core rays 
a reflecting boundary 7 + = 7 is used. Noting that v = 7 + - u, 
one gets 



Au 

d7 

for core rays and 
Au 



(p<R) 



At 



(p > R) 



(4) 



(5) 



for non-core rays. Using a suitable discretization sc heme (e.g., 
the standard scheme as modified by Pauldra ch et al.l (12001): see 
also Hoffmann et al. (2006), in prep.), the numerical solution of 
this system can be carried out very efficiently. 

A very similar method can be used to obtain a solution for 
the moments equation in spherical symmetry 



At a x ' 



77 



A 2 (qfJ) 



At 



t2 



djqfJ) 
At 



(6) 



r 1 rl 

where J = J Up) Ap, 77 = J I(p) p Ap, and all symbols with 
a tilde denote the respective quantity times r 2 . The Eddington 
factor f v is defined as 



K j\l(p)p 2 Ap 



(7) 



(at each frequency point v); q denotes the sphericality factor, 
defined by 

d(r 2 q) 



Ar 



2 3/ 
:= f-q 



1 



rf 



2 I f r 3 f{r') -\ A , 

qr =exp U ^r dr 



(8) 
(9) 



One advantage of the moments equation, as compared to 
the ray-by-ray solution, is that one can implicitly solve for the 



contribution of Thomson scattering to the source function 5. 
Separating the emissivity due to true processes [r] tme ) and that 
due to Thomson-scattering (r] Th = ^ Th 7), one can write 



x Jh J 



s = 



„Th 



r ut +x 

with the definitions 



= (1 -/3)S tme +/3J 



and 



r 



X 



Th 



X tme +X Th ' 



Using this in Eq. © gives 
A 2 (qfJ) 



At 2 



7(1 -P)-S nue (l -P). 



(10) 



(11) 



(12) 



At the boundaries, the system is closed by employing factors 
^u{p)pAp 



Jo u(p) Ap 



(13) 



similar to the second Eddington factor 1 , with u(p) coming from 
the solution of the ray-by-ray solution. At the outer boundary, 
because «(r max ) = v(r max ), this is 



h(r = r max ) = — 



Thus, the outer boundary equation is 
d(fqJ) 



At 



■ Mr : 



(14) 



(15) 



The inner boundary (r = R) is treated similarly; however, be- 
cause J upAp + 77, 7 C0le from the ray-by-ray solution has to be 
employed here as well, noting that 

ff(Tmax) = I VpAp= I I COTe pAp- I UpAp. (16) 

Jo Jo Jo 



This results in 
difqJ) 



At 



R 7 core p Ap - h(T max )J. 



(17) 



In the final, detailed solution, the radiative transfer is solved 
through the total hemisphere using a formal integral on an 
adaptive micro-grid. This also requires the specification of 
boundary conditions. For consistency, the 7^ values used for 
the sampling method must be used for the core rays in the iter- 
ations using this method. Therefore, only the boundary condi- 
tion for the sampling iteration will be discussed here. 

3. Physical conditions at the photosphere 
of a SN la 

In normal stellar atmospheres the exponential increase of the 
density at the bottom of the atmosphere provides a clear def- 
inition of a photospheric radius because large optical depths 



1 The second Eddington factor is actually defined as the ratio 
(H/J)\ r=r at the outer boundary of the atmosphere. 
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are reached at all wavelengths within a very short spatial dis- 
tance. In contrast, SN la do not have a clearly defined photo- 
sphere: because the material is unbound, no exponential density 
structure comparable to a stellar atmosphere can develop. As a 
result, the optical depth scale depends strongly on the wave- 
length. (This makes the concept of a mean optical depth like 
the Rosseland optical depth much less useful, if not entirely 
meaningless, in SN la.) In addition, the absolute densities are 
much lower than in stellar atmospheres and the composition in 
SN la is dominated by intermediate-mass and iron-group el- 
ements. When compared to stellar atmospheres this behavior 
leads to very low number densities of ions and electrons, result- 
ing in a significantly weaker free-free continuum. The absence 
of hydrogen and helium in the ejecta further reduces the contri- 
bution of the bound-free continuum in the optical and infrared 
part of the spectrum. Going from red to blue wavelengths the 
first strong continuum edge is the O i ionization edge at 91 1 A. 




1000 2000 3000 4000 
Wavelength [A] 



5000 




2000 4000 6000 _ 8000 
Wavelength [A] 



10000 



Fig. 3. Comparison of the different contributions to the total 
opacity near the photosphe re in a stellar atmosphere (upper 
panel; model D40 in Pauldr ach et all20 01 ) and a SN la model 
(lower panel); note the logarithmic scale. ^ true represents the 
total true opacity (continuum and lines), Xmnt the true contin- 
uum alone and x±om the Thomson scattering opacity. The dot- 
ted lines show the shape of a blackbody spectrum at the respec- 
tive temperature on an arbitrary scale to indicate the position of 
the flux maximum. 

Fig- E shows a comparison of the contributions to the to- 
tal opacity in the photospheric region of a hot star (Model 



D40 from Pau ldrach et al J (120011) 1 and a SN la. The dotted line 
shows a blackbody spectrum corresponding to a typical temper- 
ature (40 000 K for D40, 12 000 K for the supernova) to indicate 
the approximate position of the maximum flux in the spectrum. 
It can be seen that in the supernova, the bound-free continuum 
opacity is irrelevant compared to the line opacity in the major 
part of the spectrum. The plot also illustrates the formation of 
the "pseudo-continuum" by the overlap of thousands of lines. 

The most significant qualitative difference between a su- 
pernova and a star, however, is that at wavelengths redward of 
about 5000 A, the electron scattering opacity becomes the dom- 
inating source of opacity, even in deep layers of the ejecta. Note 
that this situation cannot be changed significantly by comput- 
ing down to smaller radii because the mild increase of density 
does not permit the formation of a significant free-free opac- 
ity. On the contrary, one encounters a lower contribution of 
true opacities because higher ionization stages tend to have less 
lines, and thus the line opacity decreases inwards. (In addition, 
the validity of the stationarity assumption needs to be consid- 
ered because the photon trapping time at deep layers becomes 
comparable to the escape time.) Fig.|4]shows the logarithm of 
the total opacity as a function of velocity and wavelength for a 
SN la model (epoch: 25 days after explosion) where the effect 
of decreasing line opacity can be seen clearly. 




Fig. 4. Logarithm of the total opacity in a SN la model (sam- 
pling iteration) versus velocity and wavelength. Note that the 
line opacity decreases toward the inside (front) because higher 
ionization stages with less lines dominate. 

Compared to a blackbody spectrum at the respective tem- 
perature, the radiation field in the innermost regions is more 
likely to have a bluer characteristic because it is the result of 
radiation from the deposition of y rays that is not entirely ther- 
malized. Furthermore, no emission from down-scattering of y- 
photons can be generated further out in the ejecta in wavelength 
regions that do not have significant continuum or line opacity. 
This effect results in the characteristic shape of SN la spectra 
in red and infrared wavelengths where the slope of the pseudo- 
continuum is generally steeper than the slope of a correspond- 



D. N. Sauer et al.: Non-LTE radiative transfer models of SN la 



7 



ing blackbody spectrum. In synthetic spectra computed assum- 
ing thermalization this effect generally results in an offset of the 
mod el spectrum in th e re d and infrared wavel engths (see e.g., 
IPauldrach et alJll996l an d iNugent et alll 19971 and the spectral 
fits in lStehle et alJ2005l) . 

In summary, the ejecta of early SN la form an intermediate 
object between an extended stellar atmosphere and a planetary 
nebula. For both extreme cases the choice for the boundary con- 
ditions is clear: for the star, the LTE diffusion approximation 
Eq. (11 8l i is a suitable choice. For a gaseous nebula, the incident 
radiation field from the illuminating star naturally defines the 
radiation field at the inner boundary (because there is essen- 
tially no back-reaction of the nebula on the stellar atmosphere). 
In SN la neither of these choices can be strictly applied. 

In the next section we will discuss an extension to the dif- 
fusion approximation that eliminates most of the restrictive re- 
quirements of LTE and therefore allows a more consistent de- 
scription of the inner boundary for supernova conditions. 



4. An improved inner boundary for the radiation 
transfer at non-LTE conditions 

4.1. The standard case for stellar atmospheres: 
the LTE diffusion approximation 

The numerical solution of the radiative transfer equation re- 
quires a set of boundary conditions. For the objects discussed 
here, this requires an assumption about the incoming radiation 
Iy at the core. 7+ should be chosen so that it describes the ra- 
diation field as accurately as possible under the physical condi- 
tions present. Ideally, the expression for the boundary equation 
is an analytic extrapolation of the radiation field at the inner- 
most points. 

For the calculations done here, the incident radiation from 
outside the ejecta is assumed to be zero. Therefore, the bound- 
ary condition for the outer boundary is /" (v, fj)\r=r m „ = 0. The 
common choice for the 7+ at the inner boundary in stellar atmo- 
sph eres is derived f rom the LTE-diffusion approximation (see, 
e.g.. lMihalasll978h 



dB v (T) 
dry 



(18) 



This result is more generally obtained by applying a har- 
monic expansion for I + that leads to the expression (see, e.g., 
IPomraninJ 19731) 



7 + (v,^) « 7°(v) + 3///V) = Iy + 3A*#? 



(19) 



where the zeroth term /J? is isotropic and the first term has 
an angular dependence proportional to ft. The term 7 1 has the 
properties of a flux (it is, however, in general not identical to 
the real flux H v that results from a solution of the transfer 
equation with this boundary condition). We therefore associate 
the symbol 77" = I (v) with the analytical expression for the 
anisotropic term of the boundary equation. In the following we 
will discuss expressions of different refinement for the terms 7° 
and 77°. 



For the classical derivation of the LTE diffusion approxima- 
tion, commonly used in stellar atmospheres, a Taylor expansion 
of S v in the limit of large t 



S v (t v )*B v = } j — ) — 



n=0 



(20) 



is used. Through the formal solution of the transfer equation 
this leads to the terms 



V - 'classic^) = B V (T) 



and 



m » h: 



classic 



(v) = 



1 dB v (T) 
3 dr 



1 1 dB y dT 
~3~x^~dT~d7' 



(21) 



(22) 



In the standard implementation, the (lower) boundary equa- 
tions for the ray-by-ray solution (core rays), Eq. ® and the 
moments equation Eq. dl7> . are therefore given by 



du v 
d(f v q v Jy) 



df v 



fy + 3(XH°y ~ Uy(T v 

(\l y+ H^)R 2 -hy(T v nJy(Ty 



(23) 



(24) 



The expansion Eq. (I20> and, therefore, Eq. ( 12 li is applica- 
ble if the radiation field is thermalized (i.e., the mean free paths 
of photons are much shorter than any significant hydrodynamic 
length scale). In SN la, however, this condition is not fulfilled 
over the full spectrum (as discussed in Section^. Hence, the 
use of Eq. (12 1 1 leads to incorrect spectral properties of the ra- 
diation field at the inner boundary. The inconsistency between 
the radiation field and the physical state of the matter caused by 
enforcing a thermal radiation field at the core boundary leads to 
spurious results in the rate equations and, in particular, in the 
heating and cooling rates for the temperature determination. 
The unstable behavior results from a feedback effect between 
the radiation field and the temperature. The radiation field at 
the inner points is used to determine the temperature, which 
(through the boundary equation) sets the incoming intensity 
and, therefore, partially sets the radiation field at those points. 

Part of this work focuses on deriving an analytical expres- 
sion for the radiation field at the inner boundary that reflects 
the physical conditions in SN la more accurately and repro- 
duces the slope of the pseudo-continuum in the red and infrared 
wavelengths better. 

4.1 .1 . Flux constraint 

To constrain the total flux at the inner boundary, the frequency 
integrated input flux H° - J H®dv is compared to the to- 
tal integrated input flux <H§ = L/(l6n 2 R 2 ). This results in a 
(frequency-independent) scaling factor 

Ho 



FC = 



ftf°dv 



(25) 



that is applied to the flux term of the inner boundary condition, 
so that (cf. Eq.fJD 



"classic W- dT d /L classic 



(26) 
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with 



FC 



classic — 



L 

\(m 2 R 2 



dT r°°j_ 

~ dr Jo Xv 



1 dB r (T) 



dv 



(27) 



(Mihalas 19781 p. 252). The purpose of applying this flux cor- 
rection factor is to induce the model to converge to a solution 
with the correct total flux. The factor FC not only accounts for 
small numerical differences (in particular in the derivation of 
the temperature gradient) but also for inconsistencies between 
the assumptions implied in the construction of the boundary 
flux and the physical conditions that are actually encountered 
at the boundary layer of the model. For the case of the clas- 
sical diffusion approximation, deviations from FC = 1 in the 
converged model indicate departures from ideal LTE diffusion 
conditions. 

With respect to the moments equation for the flux at the 
inner boundary the derived flux is actually 

H V (R) = f (/? + 3fiH?)fidfi - h v J v = H° v + (i/ v ° - hyJy) . (28) 

*J 

Thus, constraining the flux by adjusting only with respect 
to H° implicitly assumes that the second term in Eq. Il28t van- 
ishes, which requires 



(29) 



As can be seen in Fi g. [5] even for the D40 star model (cf. 
Pauldrach et al. 1200 lh this condition is not exactly fulfilled. 
In this case, however, the deviations are certainly negligible, 
whereas in SN la the effect is much stronger (bottom panel 
in Fig.[5}- This leads to a significant discrepancy between FC 
and the actually derived flux. The modifications discussed in 
the next section require that this effect is taken into account. 
Therefore, the correct flux correction factor has to be used: 



FC 



l6n 2 R 2 



(30) 



Without the additional term that accounts for small deviations 
from LTE the integrated flux at the inner boundary deviates 
from the desired input flux even if FC = 1 . With the formula- 
tion in Eq. d30i . it is possible to achieve FC = 1 and the correct 
flux at r = R. 

As already noted above, in the formulation of Eq. 1261 the 
factor FC effectively represents a correction to the temperature 
gradient at the inner boundary: 



1 8B V ( dT\ 
t (v ,,) = B v (T) + ,--(FC-). 



(31) 



We have found that it is numerically more stable to determine 
the temperature T\ of the innermost grid point directly from 
requiring that the correct flux be reached, instead of applying 
FC as above and waiting for the temperature-correction mech- 
anism of the code (based on heating and cooling rates) to ad- 
just the temperature accordingly. One reason for an unstable 
behavior of the latter approach is that in Eq. (13 1 1 the first term 
is also temperature dependent, and, in case of the physical con- 
ditions in a supernova, is of the same order as the second term. 



4x10 



S° 2x10 



DO 



-2x10" 5 



-4x10" 5 




100 



1000 



Wavelength [A] 




1000 



Wavelength [A] 



10000 



Fig. 5. The ratio {\B V — h v J v )/B Y at t he inner boundary (se e 
text). Upper panel: O-star model D40 dPauldrach et alJl200ll) . 
Lower panel: SN la model. Note the different y-scales in these 
plots. 

Therefore, a correction of the second term alone can lead to 
an inconsistency between the shape of the radiation field that 
results and the actual temperature. (An additional problem un- 
der near-LTE conditions (e.g., at the inner boundary in "nor- 
mal" stellar atmospheres) is that determining the temperature 
from the balance of the kinetic heating and cooling rates be- 
comes inaccurate, because in thermal equilibrium the heating 
and cooling rates balance at all temperatures. In stellar atmo- 
sphere models we avoid the problem by computing the temper- 
ature in the optically thick inner region either via the condition 
of radiative balance, or as a parametrized function of the op- 
tical depth, with th e param eters adjusted to conserve the flux 
(cf. IPauldrach et aljEoOll) .') From Eq. ( I30i and the condition 
FC = 1 the temperature of the innermost point is derived by 
requiring 

r^ dv -(i6^-r^-^) dv )=° (32) 

where, for the case of the standard diffusion approximation, 
corresponds to the flux term H°. . (v) of the boundary condi- 

classic v ' J 

tion (Eq. \22Y ). The conditional equation Eq. M2\ implicitly de- 
fines T\ — T2 + (dr/dr)Arj 2 because 1^ and //J? are functions of 
T and dT/dr (Eas.l21land l22t . The temperature T2 of the sec- 
ond grid point (and all other grid points) is derived convention- 
ally via heating and cooling rates. In practice, Eq. d32l > is solved 
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numerically in each iteration, assuming that the temperature- 
dependence of the second integral is negligible. 

Note that the expression Eq. d28t implicitly contains the as- 
sumption of thermalization. Only then is the degree of isotropy 
given that is necessary to make the expansion Eq. dl9> mean- 
ingful, neglecting quadratic and higher order terms. 

4.2. 1 + for a non-thermal radiation field at the 
inner boundary 

In this section we will consider modifications to the inner 
boundary that allow deviations of the radiation field from ther- 
mal equilibrium conditions, which reflects the physical situa- 
tion in SN la better. A consistent treatment of the boundary 
becomes increasingly important for models of later epochs, as 
long as the luminosity emitted at the boundary is still signif- 
icant compared to the flux originating from the y-ray energy 
deposition above that boundary. 

All modifications must be carried out in such a way that, in 
the limit of LTE-conditions at the inner points, the standard dif- 
fusive boundary condition Eq. dl 81 is retained. While it will not 
be possible to determine a boundary condition entirely free of 
analytical approximations (the model would not be sufficiently 
constrained) some of the assumptions entering into Eq. fl21i 
can be relaxed without affecting the stability of the solution. 

Starting from Eq. dl9l one can give up the assumption 
of strictly thermal conditions by allowing deviations of the 
terms / (v) and I l (v) from the Planck function. Thus, instead 
of Eq. Jl 8i we set more generally 



(33) 



with an intensity term /{' and a flux term 77° to be determined. 
We will consider those two terms separately because we will 
choose a fundamentally different treatment. 

4.2.1 . The isotropic term /{? 

We have found that the isotropic term J° can be determined 
numerically by a simple iteration procedure. This does not in- 
volve very much additional computational effort compared to 
the standard solution for the boundary condition because the 
iteration is carried out to determine the intensity term for the 
Thomson scattering source function anyway. 

In the solution of the moments equation, one can solve for 
this term implicitly by writing the boundary equation as 



dr v 1 
d(f v q v J v ) 



(34) 



df v 



(\-h v )J v =H° v R 2 (T v = T™ ax ). (35) 



An equivalent expression for the ray-by-ray solution is not as 
straightforward and may lead to slightly inconsistent results for 
J v in the iteration because there is no constraint that requires 
the angular integral of the converged analogous quantities u v 
for each ray to equal J v . For the boundary equation in the ray- 
by-ray solution we therefore also use the J v that results from 
the moments equation. 



Compared to the standard diffusion approximation this ef- 
fectively means that the inner boundary is less strongly con- 
strained, which may lead to numerical instability. To ensure 
that the boundary condition is still well behaved and to under- 
stand its general behavior we have studied this modification on 
a simple toy model before applying it to the radiative transfer 
code. 

The results of this toy model are shown in Fig.|6] The basic 
parameters are the same for all cases considered. The Thom- 
son-opacity is XTh = 1 and the background opacity (continuum) 
is Xc = 10~ 3 with a corresponding source function of S c = 1. 
We also consider a line with a source function of S\ — 5 at a 
certain radius point. (The idea here is that the line is Doppler- 
shifted by the velocity field and appears at the frequency point 
being considered at that particular radius point.) The zeroth 
term of the traditional boundary condition is set to / = 10 
for the first three cases and to 7° = 2 for the last. 77° is set 
equal to in all cases. The radiation field is obtained by an it- 
eration of a ray-by-ray solution with a solution of the moments 
equation in spherical symmetry. A Feautrier scheme similar to 
the one in the radiative transfer code is used. Generally, in the 
main code, convergence is obtained within less than 15 itera- 
tions - depending on the physical conditions and on the rel- 
evance of Thomson-scattering in particular. (For comparison, 
the iteration for Thomson-scattering alone with a fixed bound- 
ary usually converges within 5 iterations.) All plots show the 
comparison of J as a function of radius obtained from the so- 
lution of the moments equation. The result from the traditional 
choice of using a pre-specified 7° = const (in practice B v ) is 
shown as a dotted line in red. The result of allowing 7° to con- 
sistently converge to 7° = J is shown as a solid blue line. The 
black dash-dotted line represents the true source function (line 
and background). 

Panel a in Fig. |6] shows the situation for an optically thin 
model (Tii ne = 0.001). Panel b shows an intermediate case with 
Time = 0.2. The last two panels, c and d, show the case of an 
optically thick line (rii ne = 5). Here the modification only in- 
fluences the radiation field in the inner region as the emergent 
radiation is entirely separated from the inner region by the op- 
tically thick line. Compared to the conventional treatment with 
a small 7°, in the new treatment J is significantly larger in- 
ner region. The physical conditions cannot cause this increase 
in intensity because of the low true opacity. Thus, in this sit- 
uation a shortcoming of this method becomes apparent: one 
would expect J v to drop to \S\ toward the inner boundary be- 
cause the absence of emission toward the inner region means 
that 7 + = 0. The model, however, effectively sets 7 + = 7 . 
Unfortunately, this situation occurs quite frequently in SN la: at 
each frequency point where the opacity at the boundary is low 
but increases outward as a line is shifted into that frequency by 
the large velocity gradient. 

This behavior follows from the assumption that the condi- 
tions and the radiation field J v within the inner zones of the 
computational grid are representative for the region below the 
innermost point. In cases where a strong line is present further 
out this assumption is, however, not justified. If this is not taken 
into account by an additional correction, an artificial emission 
will build up in the iteration between the moments equation 
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and the formal ray-by-ray solution. Eventually this additional 
emission also affects the rate equations and the temperature in 
the inner region because the respective line transition can pump 
itself in an unphysical way. 

As a first step for an ad hoc correction of this behavior, 
criteria have to be established to determine when a correction 
should be applied. No correction is needed if the continuum 
is optically thick or if the local opacity at the first two radius 
points is large (e.g., for a strong continuum or if a line is present 
at the inner boundary). Another criterion has to include a com- 
parison of the local opacity to an average opacity over a ref- 
erence A/?-step. If the average opacity is higher than the local 
opacity, a line is likely to be present further out. The reference 
AR is chosen according to a step Ar c > 3 for pure continuum 
opacity (true and Thomson). The exact conditions for the cor- 
rection used in the current implementation are listed in Tabled 

Secondly, a suitable correction has to be used for each fre- 
quency point. We found that a suitable approach was to use a 
fixed value /J = J v > with J v > being J v of the previous (redder) 
uncorrected frequency point. To prevent excessively large val- 
ues for in frequency regions where many subsequent points 
have to be corrected, an upper cut-off at B Y is applied. 

4.2.2. The flux term H% 

Instead of the expansion Eq. d20b of S v , we now start from 
a general expression for the source function, which explicitly 
takes a contribution from Thomson scattering into account 



S v = (1 -/3 V )B V +/3J V with B y = 



Xy 



xl h +xT' 



(36) 



From the moments equation in spherical symmetry, Eq. (|6j, one 
then gets 

^ = »*(«) (37) 

(with the Eddington factor f v = J V /K Y , cf. Eq. @) which can 
be solved analytically by rewriting it to 



d 2 



Z2 (ivfv (Jy ~ By)) = 1 2 J " (<7v/v (jy ~ By)) 



df 



under the assumption that 

d 2 q v fyBy 1 ~ By 

= and — - — ss const. 



qlfy 



(39) 



The first assumption can be justified by considering only up to 
first order terms in an expansion of q Y f v By in f v . The second 

Table 1. Conditions for the correction of enhanced J v at the 
inner boundary. The reference radius R re f is set to the radius 
of tr oss ss 2/3. R Te f2 refers to the radius where t c+ i > 3. All r 
values are derived radially from the inside outward. 

Continuum: T c (/? re f) < 10 

Local*: xW(R(?)-RW) < 6 

Compare line to continuum: Tc+i(Rref2)/T c (RKo) - 3 

Compare local to average %: T c+ i/(y(l)(/? rc f - /?(!))) > 1 



assumption is not strictly fulfilled, however in practice a rep- 
resentative mean value ((1 - B v )/{qlf v )) is used. The general 
solution for J Y in Eq. J38i is then derived to 



qyfyJy = q V fyB v + C v e V \ ' + C' v 



B v + — —re 



l-A. 



a 



q v fyr 2 



qyfyr 2 



(40) 



with integration constants C v and C[, to be determined. Given 
the condition that f — > oo, J v = B v has to be obtained. It follows 
that C' v = 0. 

This result can be used to determine the flux term H® from 
the moments equation Eq. (|6j, which leads to 



H° v = —(q v f v r 2 J v ) 



= —{qyfyr 2 B v )- Cy 




For the first term one derives the expression 



(41) 
(42) 



— [qyfyr 2 By) 



d(r 2 q v ) 2 (df ^ f dB v \ 

"^/vfiv + r,^— B v + f v — j 

/ / 3/y-l d/ V \ dBy\ 2 

= \ -\ -7Z- \By + f,— \r q v 

\ \ qyXy r dr v/ dr v J 



3/v - 1 df, 



Xy r 



dB, 



(43) 



making use of the definition of the sphericality factor q v 
(Eq. (|5))- Note that the last line in Eq. J43I contains only deriva- 
tives in t v , not f v , because q v cancels in df v = -^ v * v dr. 

Next we consider the integration constant C v in Eq. (I40i . 
This constant can be obtained by considering that in the outer 
part of the atmosphere 



Jy(TQ) = — 

Jy 



for T (v) <k 1 



holds. It therefore follows that 

dqyfyJy 



(38) Jv(Tq(v)) = jy(jo) 



dfy 



(44) 



(45) 



Inserting the result of Eq. i4-0\ into Eq. ( 145 \ gives 

Cy 



Bv(ro(v)) + 



(qyfy)f 



itfy 



d [q v f v By) 




dfy 

and hence the expression for C v where to <k 1 

d(g v / v fi v ) 



(46) 



-By(f (v)) + ;'y(f (v)) 



c 



(47) 



In practice, fo has been chosen such that, at the corresponding 
depth point, the radiation field is not entirely decoupled from 
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matter (e.g., ?o ~ 0.1). This is necessary because the tempera- 
ture in B v (T(tq)) at this depth point still has to be meaningful 
to characterize the radiation field. Using the result of Eq. ( I43> . 
C v can be expressed as 



C v = 



^o) 2 {iv(g-^) 



I 



fiv(T )+/v(f ) 



T (V) 



(48) 



As a further approximation we adopt an expansion in 1/t v for 
f v assuming that f v — » | for t v — > oo to avoid mixing terms of 

f and/ v inEq.(ED: 



/v(l/Tv) 



d/v 



;(1/Tv) 



/v(Ty) * | - ^T V . (49) 



d(l/r v ) 

Introducing Eq. ( 1491 and Eq. fl43i in Eq. d41t gives a new ex 
pression for the first term of the inner boundary condition 2 




1 -ygy 
2y/v 



(50) 



This expression resembles the original diffusion flux ■ (v) 
plus additional correction terms that vanish for large optical 
depths and an isotropic radiation field. 

In the radiative transfer code, a slightly different form was 
used by solving Eq. ( 1491 1 for ^ because the gradient of the 
Eddington factor is numerically less accurate than f v itself 

(51) 

This leads to the alternative form of 




tfifv 



'-ft- 



(52) 



In the equation for the temperature at the inner boundary 
Eq. ( I32^ accordingly the expression #„ ew ,(v) has to be used 
for H®. The Eddington factor f v is calculated from the mo- 
ments J v = u v (fj) dju and K v = j Q u v (jj) /j 2 dfi of the radiation 
field, obtained from the ray-by-ray solution (which represents 
an angular-resolved computation of the intensity, I v (fi)). Note 
that f Y is a quantity not directly imposed by the boundary con- 
dition, which only specifies / + but not I~. 

The crucial modification with respect to the classical for- 
mulation Eq. (1221 of the flux term is achieved by the Eddington 
factor f y that can deviate substantially from its value of i in the 
LTE diffusion limit. This deviation can be seen in Fig.0which 
shows f v as a function of wavelength for the four innermost 
radial grid points of a SN la model. 

2 In the equation for the temperature at the inner boundary point 
(Eq. 1321 ) accordingly the expression H® cvl (v) has to be used for flj. 



5. Discussion and test calculations 

Fig. [S] shows the radiation field J v at the innermost grid point 
for the case where the new method is used with and without 
correction (indicated by the blue solid line and the green dotted 
line, respectively). The third model shown in Fig. [8] uses the 
standard boundary condition (red solid line). One can see that 
the uncorrected new boundary treatment may create large arti- 
ficial emission peaks. These peaks occur at wavelengths where 
a line is present further out within a small r-interval. 

Additionally, Fig.[8]clearly shows that the characteristic of 
the radiation field is far from Planckian, which causes the stan- 
dard diffusion approximation to be inappropriate. Also it can 
be seen that the new method produces less radiation in the red 
and infrared regions compared to the old boundary treatment. 

Looking at the structure of the newly derived flux term 
Eq. \5 2i . we note that the original flux term of the diffusion 
approximation Eq. d22t > is obtained in the limit of large r and 
f y —* h, equivalent to the requirement of a isotropic radia- 
tion field. Under these conditions the radiation field approaches 
LTE and the iterated I® = J v term will therefore approach the 
Planck function B V (T). This behavior complies with the re- 
quirement that the original diffusion approximation has to be 
recovered for LTE conditions. 

Fig. |9] shows a comparison of two full test models with 
the old and the new treatments of the boundary. One can see 
that in the red wavelength region as well as in the peaks of 
the spectrum, the radiation field in the model using the new 
method is slightly diminished compared to the model using 
the standard procedure. A direct comparison of the two mod- 
els, however, is difficult because the treatment of the boundary 
condition has a significant influence on the occupation num- 
bers, the ionization fractions, and the temperature structure. 
This influence can be seen especially in the UV part of the 
spectru m. The observed spe ctrum of SN 1992 A 5 d after max- 
imum (Kirsh ner et al.l fl993 L is shown in gray in this figure. 
Note that even though these models are test cases only and 
have not been tuned to fit the observation in detail, the overall 
shape of the observed spectrum is reproduced quite well. The 
model sh own here uses the den sity distribution of the model fl 
by Ropke & Hillebrandt (2005) which has been averaged over 
angles to obtain spherical symmetry. For the velocity field a 
homologous expansion law v oc r is assumed. For simplicity, 
we adopt a 'generic', homogeneous composition, independent 
of the predictions of the underlying explosion model. Also, the 
luminosity was not determined from the 56 Ni content of the ex- 
plosion model but set to L — 1 x 10 43 erg s . The entire lumi- 
nosity is emitted at the lower boundary implying conservation 
of the radiation flux through the ejecta. Table[2]summarizes the 
parameters of this model, vo is the velocity at the innermost 
radius of the computational grid. The velocity at the 'photo- 
sphere', v p h (where the photosphere is defined as the point at 
which r Rosse i an( ] = 2/3), is actually an output quantity of the 
model because it depends on the opacities that change with the 
occupation numbers over the course of the iterations. Due to 
the strong wavelength-dependence of the opacities, however, 
this 'photospheric velocity' is not necessarily an observation- 
ally meaningful quantity. We note that the absorption minimum 
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Table 2. Model parameters for the test model shown in Fig. [5] 
The composition is given in fractions by mass. The photo- 
spheric velocity v ph denotes the velocity where T Rosse i an d = 2/3. 
vo is the velocity of the actual innermost point of the model 
while Vph is an output of the calculation. 



C 1.1% 




Ti 


0.1% 


O 11.4% 




Cr 


0.1% 


Mg 10.2 % 




Fc 


3.1% 


Si 41.4% 




Co 


11.4% 


S 13.4 % 




Ni 


2.1 % 


Ca 5.7% 








L 1 x 10 43 


erg s _I 




Epoch 


25 


d 




v 4.48 x 10 3 


kms" 1 




Vph 6.65 x 10 3 


kms~' 





of the Si ii /16355 feature in the model corresponds to a veloc- 
ity of ~ 8 800 kms -1 which is significantly larger than the v p h 
defined via TR OSSe iand. This should be kept in mind when us- 
ing v.si n to track the photospheric velocity observationally (e.g., 
Bene tti et al.l2005tlHachinger et all2006l) . 

6. Conclusions 

We have shown that the physical conditions in the expanding 
atmospheres of SN la are such that even at early times a com- 
plete thermalization of photons cannot be assumed. Therefore 
the common approach using an LTE diffusive boundary condi- 
tion at the inner point of the computational grid does not pro- 
vide a consistent description of the radiation field and may lead 
to non-realistic synthetic spectra. Observationally this can be 
seen in the red and infrared wavelength bands where the spec- 
tral slope of typical SN la spectra deviates significantly from a 
thermal continuum. The assumption of such a thermal contin- 
uum in radiative transfer models generally results in an over- 
estimate of the radiation flux in those wavelengths. We have 
developed a theoretical framework to eliminate some of the re- 
strictions that are imposed by the assumption of LTE condi- 
tions. 

With the formalism discussed in this work we are able 
to self-consistently derive the isotropic term of the boundary 
condition. Only the flux term has to be specified analytically 
by taking the physical conditions in the photospheric region 
into account. Of course, removing constraints in the bound- 
ary conditions may cause the system to become numerically 
less stable. For this reason the consistency of the solution is 
much more important than in the classical case where the ex- 
plicit specification of the Planck-function forces the system 
into LTE-conditions naturally. 

Our modifications to the LTE-diffusion approximation have 
been derived in a very general way. Therefore, this formalism 
can also be applied to other objects (such as Wolf-Rayet stars 
with very extended atmospheres) where the physical conditions 
are such that the assumption of LTE at the photosphere is not 
justified. 



The comparison of the synthetic spectrum from a SN la test 
model to an observed SN la spectrum shows that the overall 
shape and prominent features of observed SN la are well re- 
produced by the model. A more detailed analysis of observed 
spectra will be the subject of forthcoming publications. 

Late epochs of SN la have not been considered here in more 
detail because the energy deposition by y-photons above the 
photosphere has not yet been fully implemented in our code. 
The description developed here nevertheless provides a basis 
for a reliable implementation of this energy deposition, since 
even if a fraction of the radiative energy is created above the 
photosphere the remaining radiation has to originate from be- 
low the computational grid. However, given that the photo- 
spheric conditions are such that thermalization occurs only par- 
tially, it is impossible that the energy deposited in the outer 
ejecta will be completely thermalized. Instead, excitation and 
ionization by fast electrons and y-photons above the photo- 
sphere may provide a non-thermal contribution to the spectrum. 
This might already be an issue even for epochs around maxi- 
mum, as recent three-dimensional explosion models indicate 
extensive mixing of 56 Ni into the outer layers of the ejecta. 

Furthermore, to derive a model with a luminosity that is 
consistently determined from the 56 Ni distribution of an explo- 
sion model, the time-dependent effects of photon trapping in 
earlier epochs have to be incorporated into the boundary lumi- 
nosity (lArnettll 19821 iHoflich & Khokhlovlll996t iNu gent et all 
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l =const 
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Fig. 6. Simple toy model to illustrate the new treatment of the 
zeroth term of 7 + at the inner boundary for the situation of an 
optically thin true continuum. The dotted line represents the 
run of J if the traditional setting 7° = const is assumed, the 
solid line shows the new treatment (see text). The true source 
function is indicated by the dash-dotted line. Panel a shows 
the case of an optically thin model (line and continuum) where 
Time = 0.001. Panel b shows an intermediate case (Ti me = 0.2). 
Panel c and d show the case of an optically thick line (rii ne = 5) 
where the modification only influences the inner part because 
the emergent radiation is entirely separated from the inner re- 
gion by the line. Panel d shows a situation where the method 
would overestimate J in the inner region: the case where in the 
conventional treatment 7° at the inner boundary is smaller than 
S ii ne . In this case the iteratively determined 7° is not appropri- 
ate, and instead an 7° based on that of a neighboring frequency 
point is used. 
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Fig. 7. The Eddington factor f v = J v /K v at the innermost ra- 
dius points as a function of wavelength for the case of a SN la 
model. The deviation of f v from its LTE value of A (indicated 
by the dotted line) allows the new flux term of the inner bound- 
ary 77 n ) ew (v) to deviate substantially from the LTE value of 
77° assic (v). The radii are given in units of the innermost radius. 
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Fig. 8. The radiation field J v at the inner boundary for differ- 
ent treatments of 7 + . The dotted green line shows the radiation 
field for the uncorrected iterative method, which produces large 
artificial peaks at wavelength points where a line is present fur- 
ther out within a small r-interval. The blue solid line represents 
the new method including the correction for those wavelength 
points. The red line is calculated with the traditional treatment 
of the inner boundary. It is clear that the new method signifi- 
cantly reduces the radiation field in the red part of the spectrum. 
All three models are shown after a few iterations before the first 
temperature update. 
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Fig. 9. Comparison of two test models using the old and the 
new treatments of the inner boundary. One can see that the 
flux in the red wavelengths is diminished in the model using 
the new boundary treatment. Note that the two models are not 
strictly comparable because the occupation numbers, ioniza- 
tion, and temperature structure adjust differently. For compari- 
son the o bserved spectrum o f SN 1992 A 5 d after maximum is 
shown dKirshneret al.ll993l) . 



